Ionic solutions play a significant role in many areas including chemical engineering^[1–3] and biological systems.^[4,5] The microscopic structure and transport properties of ions are two classical topics in chemical physics. The Debye–Hückel theory^[6] describes ions as fully decomposed and solvated in solutions. However, recent studies suggested that ions tend to form clusters even in diluted solution.^[7–10] The microscopic structures of ions are essential for many processes, including nucleation,^[11,12] battery performance,^[13–15] and biological systems.^[16–18] Moreover, ion transport properties usually dominantly determine the performance of ionic solutions working in a non-equilibrium processes, e.g., driven by an external electric field E.

The transport properties, including diffusion and viscosity, are significantly influenced by the microscopic structure of solvation shells. Many works have been carried out to explore transport properties of ion under E. Lee and Rasaiah^[19] systematically studied the extremely dilute solutions by molecular dynamics simulations. The relation between solvation structure and viscosity was also explored in extremely dilute solutions with real or virtual particles.^[20] The frictional coefficients of ion were investigated under the ambient condition.^[21] It was found that the charge transport and ion-ion correlation are different for different ionic solutions.^[22] In addition, the concentrated NaCl solution under E was studied by several groups.^[23,24] Kerisit et al. exploited the dynamic and transport properties of alkali ions in dimethyl sulfoxide solutions under E.^[25]

The diffusion and viscosity are closely connected by the Stokes–Einstein relation D∼ T/η in the linear response region. However, the Stokes–Einstein relation is invalid when external applied force is large, such as under a strong E.^[24,25] Ionic solution is in a non-equilibrium steady state for the steady flux caused by E. The glass-forming liquids under undercooled state are also out of equilibrium for the slow dynamics. The Stokes–Einstein relation is observed to break down in glass-forming liquids and liquid metal.^{[13,14,26–32]} Because of the difficulty to accurately determine shear viscosity η via simulation, the structural relaxation time τ is usually adopted as a substitute to evaluate the shear viscosity,^[27,33–36] and the corresponding variant of the Stokes–Einstein relation is D ∼ τ⁻¹.^[34,35] As glass-forming liquids undergo supercooling, D ∼ τ⁻¹ is breakdown and follows a fractional form as D ∼ τ^ξ with ξ ≠ −1. Since both glassy-forming liquids under undercooled state and ionic solution under strong E are out of equilibrium, but different forms of Stokes–Einstein relation are adopted in the two cases, an interesting question arises: is there an analogous fractional variant of Stokes–Einstein relation D ∼ τ^ξ or dynamic heterogeneity existing in ionic solution under strong E?

In this work, we take three aqueous alkali chloride solutions (NaCl, KCl, and RbCl) to explore the variant of the Stokes–Einstein relation in ionic solution under strong E. The paper is organized as follows: Section 2 contains a brief description of the simulation details. The simulation results are described in Section 3. Section 4 includes our conclusions and discussion.

To explore the variant of the Stokes–Einstein relation in ionic solution under E, we choose three aqueous alkali chloride (NaCl, KCl and RbCl) solutions. Each system contains 128 ion pairs and 3072 water molecules, the concentration is ∼ 2.2 mol/L. Fourteen different E’s uniformly distributed among 0–1.3 V/nm are applied in the X direction. The diffusion and viscosity of alkali ions have been studied by Lee et al.^[20] In this work, we choose the same model for alkali ions and SPC/E water model. The related force field data are listed in Table 1. All our simulations are carried out with the GROMACS package,^[37,38] in which the external electric field took effect on a charged atom or ion (with a partial charge q) by applying an additional electric force F_e = qE in the X direction besides the routine MD forces in equilibrium, and the equation of motion for atom or ion i with mass m_i and charge q_i is

Table 1.

Table 1.

Table 1. The force field parameters for the Na+, K+, Rb+ and SPC/E water model. . Atom/ion σio/nm εio/(kJ/mol) Charge/e Na+ 0.2876 0.5216 +1 K+ 0.325 0.5216 +1 Rb+ 0.3348 0.5216 +1 O 0.3169 0.6502 −0.8476 H – – +0.4238

Table 1. The force field parameters for the Na+, K+, Rb+ and SPC/E water model. .

Ions exists as solvated ions in aqueous ionic solution. To characterize the influence of E on solvation structure, we adopt the radial distribution function (RDF).^[19] The RDF is defined as g(r) = 〈 δ (r–r₁) 〉, where r₁ is the position of the central ion. The RDFs for Na⁺–Cl⁻, K⁺-Cl⁻, and Rb⁺–Cl⁻ solution show almost the same changes under E. We only plot the RDFs for Na⁺–Cl⁻ solution to illustrate the influence of E and the related results are shown in Fig. 1.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Radial distribution functions (RDFs) of solvation shells under different E: (a) Na+–Cl− shell, (b) Na+–H2O shell, (c) Cl−–H2O shell.

Figure 1(a) shows that the main peak of the Na⁺–Cl⁻ shell has a large increase with increasing E, whereas the second peak get a small decrease. The main peak is corresponding to the ion atmosphere, which mainly determines the dynamics of ions. The results suggest that Na⁺ and Cl⁻ get more strong spatial correlation in the ion atmosphere. It is resulted from the decrement of dielectric constant as the dipole of water tends to get a same orientation under E.^[42] The main peaks of the Na⁺–H₂O shell and the Cl⁻–H₂O shell have a little decrease with increasing E as shown in Figs. 1(b) and 1(c). Ion and water molecules get less spatial-correlated under a stronger E.

The variant of the Stokes–Einstein relation D ∼ τ⁻¹ is based on the exponential relaxation of the self-intermediate scattering function F_s(k,t) ∼ e^−t/τ in simple liquids, which can be described by F_s(k, t) = exp (−k²Dt) if the displacements of particle follow Gaussian distribution.^[28] To evaluate the variants D ∼ τ⁻¹, the diffusion constant of ions is calculated via the mean square displacement as

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. The drift velocities v for Na+, K+, Rb+, and their corresponding Cl− as a function of E: (a) v for Na+, K+, Rb+, (b) v for alkali ion corresponding Cl−. The v of alkali ion corresponding Cl− in (b) is plotted in the same color as shown in (a).

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The drift velocities v for Na+ and its corresponding Cl− versus E. The symbols are the simulated data. The solid lines are fitted by v ∼ E and the dotted lines are fitted by v ∼ E + a ⋅ E2.

The mean square displacements (MSDs) of Na⁺, K⁺, Rb⁺, and their corresponding Cl⁻ are calculated. To illustrate the MSDs under E, the results of NaCl solution is plotted in Fig. 4. The variation of MSDs for Na⁺ and Cl⁻ with time can be separated to three regimes similar as in undercooled liquids.^[28] The ballistic regime is at short times (∼ 0.1 ps) so that MSD is a quadratic function in time. The diffusive regime is at long times (> 1 ps) so that MSD is a linear function in time. The cage regime is between the ballistic and the diffusive regime (0.1–1 ps), and the logarithm of MSD increases nonlinearly with logarithm of time. The MSD in the cage regime is not a plateau like undercooled liquids. The ion is not so strongly trapped in the cage as in undercooled liquids. The cage regime shortens with increasing E, and the rate of increase of MSD with time is greater under a stronger E for both Na⁺ and Cl⁻. This is due to fact that the unidirectional drift velocity caused ion becomes more easily to get out of the trap of solvation shell. The diffusion constants of Na⁺, K⁺, Rb⁺ and their corresponding Cl⁻ calculated by Eq. (3) are plotted in Fig. 5. The diffusion constant D for each ion increases with increasing E. The tendency is almost the same for all ions and the only difference is the rate of increase.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. The logarithm of mean square displacements for Na+ and Cl− as a function of the logarithm of time t under different E’s: (a) Na+, (b) Cl−.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. The diffusion constants D of Na+, K+, Rb+ and their corresponding Cl− as a function of E. (a) D for Na+, K+, Rb+, (b) D for alkali ion corresponding Cl−. The D of alkali ion corresponding to Cl− in (b) is plotted in the same color as shown in (a).

The structural relaxation of ion is described by the self-intermediate scattering function, namely,

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. The structure factor S(k) for each ion under E = 0. S(k) for Na+, K+, and Rb+ are plotted in (a) and their corresponding Cl− are plotted in (b).

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. The self-intermediate scattering function Fs(k, t) for Na+ and Cl− as a function of time t under different E’s: (a) Na+, (b) Cl−.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. The structural relaxation time τ of Na+, K+, Rb+ and their corresponding Cl− as a function of E. (a) τ for Na+, K+, Rb+, (b) τ for alkali ion corresponding Cl−. The data of alkali ion corresponding Cl− in (b) is plotted in the same color as shown in (a).

The simulation can be well fitted by log (D),log(τ⁻¹) ∼ a + bE + cE². The fitted a, b, and c are listed in Table 2. Similar to the variation of velocity with E plotted in Fig. 3, the non-linear effect also plays a role in the changes of diffusion and relation. For small E, the linear item plays a more important role E in relaxation, but the quadratic item starts to play a role with increasing E. The linear item plays a more important role than the quadratic item in diffusion for all E considered.

Figures 5 and 8 show that the diffusion constant or relaxation times has similar changes with increasing E for all ions. To quantify the variation, the logarithms of the diffusion constant and relaxation time for Na⁺, Cl⁻ as a function of E are plotted in Fig. 9.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. The logarithms of the diffusion constant and relaxation time versus E for Na+ and its corresponding Cl−: (a) log(D) vs. E, (b) log(τ−1) vs. E. The symbols are the simulated data and the solid lines are fitted by log (D), log(τ−1) ∼ a + bE + cE2.

To examine the existence of the fractional variant of the Stokes–Einstein relation in ionic solution under E, the variant D ∼ τ⁻¹ is evaluated by D ∼ τ^ξ. It will be valid if ξ = −1, otherwise it will follow a fractional form. The logarithms of τ⁻¹ and D for Na⁺, K⁺, Rb⁺ and their corresponding Cl⁻ are plotted in Fig. 10. The fitted ξ are not equal to −1 for all ions under E. The data for all ions can be well fitted by D ∼ τ^ξ with different ξ, and a crossover is observed for all ions with increasing E. The crossover E_x is around 0.7 V/nm, 0.5 V/nm, 0.5 V/nm for Na⁺, K⁺and Rb⁺, respectively, E_x for their corresponding Cl⁻ is around 0.7 V/nm, 0.6 V/nm, 0.5 V/nm. The similar phenomenon is also observed in water under supercooling.^[27] The exponents ξ for each ion are all different from 1 and in a fraction, which show that the fractional variant of the Stokes–Einstein relation exists in ionic solution under E.

Fig. 10.

Figure Option ViewDownloadNew Window

Fig. 10. Test of the variant of the Stokes–Einstein relation described by D ∼ τ−1 for Na+, K+, Rb+ and their corresponding Cl−: (a) Na+, K+and Rb+, (b) data for alkali ion corresponding Cl−. The symbols are the simulated data and the solid lines are fitted by D ∼ τ−ξ. The colored exponent ξ is corresponding to the same colored solid line. The data of alkali ion corresponding to Cl− in (b) is plotted in the same color as shown in (a).

Table 2.

Table 2.

Table 2. Fitting data of the logarithm of the diffusion constant D and structural relaxation time τ by log(D), log(τ−1) ∼ a + bE + cE2. . DNa+ DCl− τNa+ τCl− a −7.020 −6.878 −0.448 −1.328 b 1.702 2.180 0.144 0.151 c −0.119 −0.358 0.237 0.341

Table 2. Fitting data of the logarithm of the diffusion constant D and structural relaxation time τ by log(D), log(τ−1) ∼ a + bE + cE2. .

The variant of the Stokes–Einstein relation D ∼ τ⁻¹ will be a strict result if the displacements of particle follow Gaussian distribution. It is proposed that the fractional form is resulted from dynamic heterogeneity in glass-forming liquids.^[47] The dynamic heterogeneity can be characterized by the deviations from Gaussian distribution in particle displacement. We calculate the self van Hove function G_s (r, t) and its first-order approximation G₀ (r, t) to describe the deviation.^[48] G_s (r, t) is the Fourier transformation of F_s(k,t) and G₀(r, t) takes the Gaussian form. If the distribution of ion displacement takes a Gaussian form, G₀ (r, t) and G_s (r, t) are the same, otherwise G_s (r, t) deviates from G₀ (r, t). The calculated results for Na⁺ and Cl⁻ in NaCl solution at t* = 10 ps are shown in Fig. 11. G_s(r, t*) deviates from the G₀(r, t*) under E. The results indicate that the dynamic heterogeneity appears under E.

Fig. 11.

	Figure Option ViewDownloadNew Window
	Fig. 11. Self van Hove function Gs (r, t*) and its first-order approximation G0(r, t*) for Na+ and Cl− in NaCl solution at t* = 10 ps under different E’s. Gs(r, t*) and G0(r, t*) are plotted in the same color with the dotted lines for Gs(r, t*) and the solid lines for G0(r, t*).

To explore the reasons for the crossover and the fractional form observed in D ∼ τ^−ξ, as data plotted in Fig. 9 and listed in Table 2, the logarithms of the diffusion constant and structural relaxation time vary with E as log (D),log(τ−1) ∼ a + bE + cE². The a item in diffusion is much larger than the a item in relaxation. The b item in diffusion is much greater than the c item for all E considered. Thus the response of the diffusion to E is mainly determined by the linear term in E. However, the b and c items are in the same magnitude in relaxation for all E considered. For small E, the b item is greater than the c item in relaxation but it is reverse at a large E. If we make bE = cE² in the relaxation, we can calculate the crossover E_x: E_x = 0.607 for Na⁺, E_x = 0.450 for Cl⁻. E_x is roughly equal to the observed linear–nonlinear crossover E_x in velocity and the crossover E_x observed in D ∼ τ^−ξ. Therefore, the fractional D ∼ τ^−ξ and the crossover are mainly due to the linear and nonlinear effect caused by E: at small E, the linear term plays an important role and the quadratic term play a more important role at large E.

The atomistic MD simulations for aqueous NaCl, KCl and RbCl solutions have been carried out to exploit the fractional variant of the Stokes–Einstein relation in ionic solutions under strong E. Alkali ion and Cl⁻ get more spatially strongly correlated under a stronger E but the hydration shell is weaker. Ion gets a larger velocity with a larger E that leads the system further away from the equilibrium state. E accelerates the diffusion and structural relaxation of ion, and ion gets a greater diffusion constant and smaller relaxation time under a stronger E. The variant of the Stokes–Einstein relation is found to be breakdown under E, a fractional form and a crossover are observed as E increases. The fractional form and crossover are resulted from nonlinear response introduced by E. The fractional form and crossover are caused by the dynamic heterogeneity of ion and the linear–nonlinear response under E. Our results indicate the similarities between ionic solutions under E and glass-forming liquids, which may help us to improve our understanding and applications of ionic solutions. Moreover, the fractional variant of the Stokes–Einstein relation is caused by E in our work. Actually, some evidences have shown that both the dynamics and structures of glass-forming liquids are heterogeneous. Thus our results may make some hints on how to explain the fractional variant and crossover in glass-forming liquids under cooling, such as the heterogeneity cause molecular mean field.